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ABSTRACT 


The  first  particle-in-cell  (PIC)  calculation  representing 
the  vortex  by  the  particles  was,  to  the  authors'  knowledge, 
done  by  Rosenhead  in  1930  in  his  study  of  vortex  streets.  The 
same  calculation  was  repeated  by  Birkhoff  and  Fisher  in  1959 
and  by  Hama  and  Burke  in  1960.  Abernathy  and  Kronauer  made 
a detailed  study  of  the  Karman  vortex  street  in  1962  using 
the  same  method.  Hockney  (1962)  was  the  first  one  to  incor- 
porate the  advantages  of  the  fast  Fourier  transform  (FFT)  to 
the  PIC  method  in  order  to  handle  very  large  numbers  of  par- 
ticles, usually  in  the  order  of  10^  in  the  galaxy  simulations. 
Chorin  (1973)  proposed  a new  version  of  the  discrete  vortex 
method  where  vortices  are  needed  only  to  satisfy  the  no-slip 
boundary  conditions.  Along  these  lines,  due  to  the  natural 
stratification  of  the  atmosphere  and  the  existence  of  wind 
shear,  almost  all  atmospheric  flows  are  rotational  as  a re- 
sult of  the  interaction  between  shear  and  buoyancy.  The  nu- 
merical study  presented  herein  combines  the  concepts  of  the 
PIC  method,  Green's  function  and  the  FFT  methods  and  may  be 
dubbed  the  vortex-in-cell  (VIC)  method. 

The  discrete  vortices  are  treated  as  the  marker  parti- 
cles on  an  Eulerian  grid  with  the  velocity  field  solved  from 
the  updated  vorticity  distribution.  The  velocity  field  can 
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either  be  calculated  from  the  distribution  of  the  vortices 
by  Green's  function  formalism  or  by  the  FFT.  The  computation 
time  required  for  Green's  function  formalism  is  proportional 
to  the  square  of  the  number  of  the  vortices  N2,  while  that  of 
the  FFT  is  fixed  to  MxMy  An(MxMy)  , where  Mx,My  are  the  num- 
ber of  mesh  points  in  the  x and  y directions.  Therefore, 
the  FFT  formalism  can  deal  with  a problem  of  a very  large 
number  of  particles  without  increasing  the  computation  time. 
The  VIC  method  can  also  simulate  the  physical  process  ex- 
actly through  which  the  flow  containing  or  generating  vor- 
ticity  evolves.  If  it  is  clear  that  the  essence  of  turbu- 
lence is  vortex  interactions  and  decays,  then  certainly  a 
method  directly  dealing  with  such  motion  should  be  a valid 
simulation  of  the  turbulence.  The  inclusion  of  the  viscos- 
ity using  Chorin's  (1973)  method,  and  generalization  to  three-- 
dimensional formalism,  will  further  extend  the  range  ci:  the 
relevance  to  turbulence.  The  basic  economy  saved  by  the  VIC 
method  is  not  only  due  to  the  efficient  numerical  aspects  of 
the  schemes  but  also  is  due  to  the  basic  idea  that  th  i flow 
is  dominated  by  the  rotational  motion  generated  by  the  dis- 
cret^zec^  vortices,  and  it  is  wasteful  to  compute  the  passive 
portion  of  the  flow  which  stays  passive  and  irrotational . 
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Various  problems  in  different  applications  have  been 
solved  by  this  method;  these  are: 


Meteorology 

Oceanography 


Aerodynamics 


Flow  in  porous 
medium 


Kelvin-Helmhol tz  wave  generation 

Collapsing  wakes  at  ocean  thermocline 

Buoyant  wakes  near  an  ocean  thermocline 

Karman  vortex  street 

Aircraft,  trailing  vortex  in  a wind 
shear  field 

A thermal 

An  injection  cylinder 

A study  of  interfacial  finger-like 
flow  structures. 


Good  agreement  with  experiments  was  obtained. 
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DEVELOPMENT  OF  VORTEX  STRUCTURES  IN  BUOYANT  AND  SHEAR  FLOWS 
I . INTRODUCTION 

The  formation  of  ring  or  line  pair  vortices  following 
injection  of  a blob  of  fluid  into  a medium  at  rest  is  a very 
common  phenomenon.  Some  well-known  examples  include  the 
classic  "smoke  ring",  in  which  a blob  of  air  is  injected  into 
a quiescent  environment  at  sane  initial  velocity.  The  two- 
dimensional  analog  of  this  axisymmetric  ring  vortex  formation 
occurs  in  the  wakes  behind  lifting  aircraft.  The  airstream 
behind  the  wings  of  the  aircraft  receives  a downward  momentum 
per  unit  length  along  the  aircraft  track  equal  to  the  aircraft 
lift.  Since  the  downward  velocity  of  the  air  is  usually  much 
smaller  than  the  forward  speed  of  the  aircraft,  this  phenomenon 
can  be  thought  of  approximately  as  equivalent  to  the  two 
dimensional  motion  that  results  when  a cylindrical  column  of 
air  is  given  an  initial  downward  velocity.  This  motion  is 
often  observed  to  result  in  two  well-defined  line  vortices 
trailing  behind  the  aircraft  wing  tips. 


When  a spherical  volume  of  air  in  the  atmosphere  is  heated 
at  constant  pressure  as  a result  of  a point  release  of  energy 
(explosion),  the  hot  bubble  of  gas  rises  and  is  typically 
observed  to  develop  into  a toroidal  ring  vortex  configuration 
by  the  time  it  has  risen  a distance  of  the  order  of  its  initial 
diameter.  The  two-dimens ional  analog  of  this  axisymmetric 
motion  can  occasionally  be  observed  in  bent-over  chimney  plumes 
where  the  hot  gas  emanating  from  the  chimney  is  stretched 
horizontally  into  a long  cylinder  by  an  ambient  wind  and  then 
rises  buoyantly.  When  the  air  is  calm  this  long  cylinder  of 
hot  gas  can  be  observed  to  develop  into  two  parallel  line 
vortices  as  it  floats  upwards. 
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In  the  case  of  the  smoke  ring  ejection  and  the  aircraft 
wake  the  initial  motion  consists  of  a distributed  sheet  of 
vorticity  embedded  in  a hydrodynamic ally  irrotational  flow. 

The  subsequent  development  of  the  motion  in  an  incompressible 
fluid  (and  in  an  ordinary  compressible  fluid  when  the  velocities 
are  small  compared  with  the  speed  ©if  sound)  can  be  thought  of 
as  a mutual  induction  or  interaction  between  the  various  sections 
of  the  vortex  sheet.  In  the  two-dimensional  case  in  the  absence 
of  viscosity,  vorticity  is  a transferable  quantity  and  no  new 
vorticity  is  created  during  the  ensuing  motion.  In  the  case 
of  the  buoyant  motions  (rising  bubble,  bent-over  chimney  plume) 
the  initial  state  may  be  approximated  as  an  irrotational  one 
with  zero  vorticity.  The  interaction  between  the  gravitational 
pressure  gradient  and  the  density  gradients  provide  a source 
of  vorticity  and  the  motion  consists  essentially  of  two 
processes  - the  creation  of  a vortex  sheet  or  a vortex  layer 
as  a result  of  buoyancy  and  the  subsequent  mutual  induction 
of  different  portions  of  this  sheet.  The  latter  induction 
phenomenon  closely  parallels  the  development  in  non— buoyant 
injection  flows.  In  the  buoyant  phenomena  the  mechanics  of 
the  vorticity  generation  are  essentially  those  that  give  rise 
to  the  well-known  Rayleigh-Taylor  instability.  The  motion  of 
the  vortex  sheet  in  the  ejection  class  of  flows  is  closely 
related  to  the  Kelv in-Helmholtz  instability  in  shearing  flows. 

A related  phenomenon  is  the  Kelvin-Helmholtz  instability 
that  occurs  at  the  interface  between  two  fluids  which  are  in 
relative  motion.  In  this  case  the  interface  (vortex  sheet) 
is  unstable  to  the  development  of  line  vortices  whose  axes 
are  perpendicular  to  the  velocity  difference  vector.  In  the 
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case  of  a wake-like  flow,  which  can  be  roughly  considered  as 
two  parallel  vortex  sheets  of  layers  of  opposite  sign,  this 
instability  results  in  the  development  of  the  well-known 
Kerman  vortex  street. 

When  horizontal  shearing  exists  in  a gravitational  field 
in  which  the  fluid  is  stably  density  stratified  (density  increasing 
upwards)  both  the  Kelvin-Helmholtz  mechanism  and  the  Rayleigh- 
Taylor  mechanism  are  operative.  In  the  case  of  stable  stratif- 
ication, the  Rayleigh-Taylor  mechanism  is  stabilizing  whereas 
the  Kelvin-Helmholtz  mechanism  is  destabilizing.  Under  these 
conditions,  instabilities  will  develop  when  the  destabilizing 
forces  dominate  over  the  stabilizing  forces.  The  character  of 
the  motion  depends  on  the  value  of  the  Richardson's  number 

g (dpT/dz) 

Rl  ~ oT(du/dz)2 

and  instability  occurs  when  Ri  decreases  below  a certain  value. 

Here  PT  is  the  density  for  incompressible  fluids  or  the  "total”  * 
density  for  compressible  fluids. 

This  phenomenon  is  presently  considered  to  be  one  of  the 
major  mechanisms  creating  turbulence  in  clear  air  at  tropopause 
altitudes.  In  this  case  the  stabilizing  influence  of  an 
upwardly  increasing  potential  temperature  in  the  atmosphere 
inhibits  development  of  the  instabilities  until  a sufficiently 
high  shearing  motion  can  develop.  The  effect  is  of  particular 
interest  since,  when  the  shear  does  develop  to  a high  enough 
level  to  destabilize,  there  is  potentially  a large  amount  of 
kinetic  energy  which  can  be  rapidly  converted  into  vortex 
formation  and  subsequently  into  turbulent  motion. 

The  conditions  for  onset  of  these  various  instabilities, 

and  characteristic  wave  numbers  and  growth  times  involved,  have 
*The  word  "total"  refers  to  the  density  of  a fluid  element 
when  brought  isentropically  to  a fixed  reference  pressure. 
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been  studied  extensively  theoretically.  Experimental  obser- 
vations are  available  in  a variety  of  laboratory  simulations 
and  in  various  natural  situations.  In  addition  to  the  examples 
previously  mentioned,  phenomena  which  are  associated  with  the 
unstable  development  of  vortex  structures  include  a wide 
variety  of  musical  wind  instruments  in  which  the  sound  is 
dependent  on  the  formation  of  vortex  streets,  the  flapping  of 
flags  in  a breeze,  the  whistling  of  wind  through  wires  and 
trees  (Aeolian  tones),  aerodynamically  induced  instabilities 
aft  of  bluff  structures  (of  which  a classic  example  is  the 
destruction  of  the  Tacoma  Narrows  suspension  bridge  in  1939), 
tornados,  dust  devils,  thunderheadu , atmospheric  thermals,  and 
hurricanes . 

There  is  another  class  of  flows  which,  though  not  always 
resulting  in  vortex  formation,  is  quite  similar  in  its  hydro- 
dynamics to  the  flows  discussed  above.  These  flows  obtain  in 
porous  media  where  the  inertia  forces  are  negligible  compared 
with  the  viscous  forces.  Here  the  mechanisms  driving  the  motion 
are  again  either  the  initial  vorticity  or  the  vorticity  induced 
as  a result  of  buoyant  or  pressure  gradient  forces.  Examples  of 
this  class  of  flows  include  the  motion  of  a variable  density 
incompressible  fluid  through  a porous  medium  under  the  influence 
of  gravity  or  a pressure  gradient.  In  the  case  of  the  gravit- 
ational flows  when  a heavy  liquid  overlies  a light  liquid,  the 
Rayleigh-Taylor  mechanism  is  operative  and  the  fluid  will  develop 
a motion  to  allow  the  heavy  liquid  to  interpenetrate  the  lighter 
liquid  in  order  to  fall  downward.  This  motion  has  an  interesting 
analog  in  certain  magnetized  low- density  plasmas.  When  a low- 
density  plasma  is  placed  in  a strong  magnetic  field  and  is 
acted  on  either  by  gravity  or  weakly  couples  to  a moving  neutral 
background  wind  (the  viscous  interaction  between  the  neutral 
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wind  and  the  ionized  plasma  generally  provides  the  coupling) . 

The  plasma  coupling  to  the  neutral  background,  however,  is 
inhibited  by  the  magnetic  field.  When  the  plasma  conductivity 
perpendicular  to  the  magnetic  field  varies  with  position, 
motion  can  develop  which  is  described  (in  the  limiting  case  of 
a very  strong  magnetic  field  and  where  the  gradients  parallel 
to  the  magnetic  field  are  negligible)  by  equations  that  are 
essentially  identical  to  the  equations  describing  the  motion 
of  a variable  density  incompressible  fluid  in  a porous  medium 
acted  on  by  gravity  or  pressure  gradient  forces.  Instability 
can  develop  to  allow  the  higher  conductivity  plasma  to  couple 
better  with  the  background  neutrals  and  allow  it  to  inter- 
penetrate the  lower  density  plasma.  This  phenomenon  is 
expected  to  be  operative  in  the  normal  ionosphere  as  a result 
both  of  ambient  electric  fields  and  winds  in  the  background 
neutral  atmosphere.  These  neutral  winds  are  collisionally 
coupled  to  the  magnetized  plasma  and  a dynamo  action  results 
which  creates  the  electric  fields  that  move  the  plasma.  It 
is  thought  that  some  of  the  irregular  structure  of  the  electron 
density  distribution  in  the  ionosphere  (which  is  manifested  as 
spread— F ) results  from  these  instabilities  (variously  called 
the  gradient  drift  instability,  the  E x B instability,  or  the 
Simon  instability) . 

A phenomenon  first  observed  in  pumping  of  oil  wells  is 
known  as  "water  tonguing"  or  "water  coning".  Here  when  the 
pressure  at  the  drill  stem  inlet  is  made  too  low  (i.e.  by 
pumping  too  hard)  water  is  observed  to  be  mixed  with  the  pumped 
oil.  In  these  cases  the  oil  is  distributed  through  a sandy 
medium  and  the  mechanics  of  the  pumping  consists  of  lowering 
the  pressure  at  the  well  inlet  so  that  the  water  surrounding 
the  oil  bed  can,  under  its  own  high  pressure,  force  the  oil 
to  the  well  inlet.  G.I.  Taylor  studied  this  effect  and  concluded 
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that  when  a low  viscosity  fluid  pushes  a higher  viscosity 
fluid  through  a porous  medium  the  interface  can  become  unstable 
at  sufficiently  high  pressure  gradients  and  tongues  of  water 
can  snake  through  the  oil  towards  the  low  pressure  point.  The 
equations  of  motion  and  the  resultant  instabilities  are  closely 
related  to  the  mechanics  of  a two-density  fluid  in  a porous 
medium,  moving  under  the  influence  of  gravity. 


In  this  note  we  consider  a class  of  flows  in  the 
special  case  where  the  flow  consists  of  two  incompressible 


' 
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II.  FORMULATION  OF  THE  PROBLEM 


II.  1 General 


The  momentum  equation  for  an  incompressible  fluid  may 
be  written  in  the  form 


du 

dt 


- — <7p  + r|V  u + g 


(1) 


— * 

where  u is  the  velocity  p the  pressure  and  p the  fluid 
density,  n the  kinematic  viscosity,  and  g is  the  acceleration 
due  to  gravity.  The  vorticity 


= 7 x u 


(2) 


satisfies  the  equation  obtained  by  taking  the  curl  of  Equation  (1) 


dt 


- ?—  x 7p  + ny-C 
o 


(3) 


Thus  vorticity  is  generated  as  a result  of  buoyancy  forces 
associated  with  density  and  pressure  gradients  and  diffusively 
dissipates  as  a result  of  viscosity.  Equation  (3)  provides  a 
means  for  evaluating  the  vorticity  of  given  fluid  elements. 

The  velocity  at  any  point  (r)  in  the  fluid  may  be  evaluated 
from  the  kinematic  identity  (in  two  dimensions) 


u (r)  = 


2tt 


// 


C(r')  x (r  - r ')  dx  'dy ' 

I V _ V ' I 2 


(4) 


where  the  integral  extends  over  the  entire  rotational  region 
of  the  fluid.  The  continuity  equation  for  an  incompressible 
fluid: 


do 

dt 


_ o_c 

fit 


+ U.vo  = 0 


(5) 


can  be  used  to  follow  the  evolution  of  the  density  distribution 
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ir  time. 


Instead  of  the  Green's  function  form  for  the  velocity 
field  (Eq.4),  the  velocity  may  be  expressed  in  terms  of  a 
stream  function  Y: 


defined  such  that  V.Y  = 0 (the  vanishing  of  v.Y  is  automatically 
fulfilled  in  two  dimensional  motion) . The  stream  function 

satisfies  a Poisson  equation  with  the  vorticity  as  the  source 
function. 


VaY  = C 


(7) 


for  which  the  formal  solution  can  be  written 


Y(r) 


In  |r-r'[  dx*dy' 


(8) 


Thf!  cur*  of  this  expression  yields  the  identity  in  Eq.  (4)  . 

In  the  present  analysis  we  will  be  concerned  primarily 

with  fluids  where  the  fractional  variation  of  the  density  and 

viscosities  are  small  (A£_,  Arp  <<;  and  we  carry  out 

P T) 

calculations  only  to  first  order  in  these  variations,  in 
Section  II. 2 we  consider  the  two  dimensional  motion  of  inviscid 
fluids  and  in  Section  II. 3 the  motion  of  an  incompressible 
fluid  in  a porous  medium,  or  equivalently,  the  motion  of  an 
incompressible  fluid  confined  to  move  between  two  closely 
spaced  vertical  walls  (Hele-Shaw  cell) . The  case  of  large 
density  differences  is  treated  in  Section  II. 4. 


I I . 2 Inviscid  Fluids 

When  the  viscosity  is  negligible  the  vorticity  equation  (3) 
takes  the  form 


K = 

dt 


1 

-v—  x 

P 


vp 


+ g 


(9) 


and  to,  first  order  in  the  density  difference  (since  p = 

-*  Ad 
o g + oef), 
o P 


d C Vo 

dt  = " T X g • (i * * * * * * * * 10> 

o 

For  two  dimensional  motion  in  the  (x,y)  plane  the 
vorticity  is  effectively  a scalar  (i.e.  has  only  a z component). 
Thus,  for  a vertical  (y  direction)  downward  gravitational 
acceleration  g,  Equation  (10)  becomes 


i = a i“  (id 

dt  p dx 

o 

Of  particular  interest  is  the  case  of  two  uniform  immiscible 

fluids  of  slightly  differing  density.  In  this  case  vorticity 

is  generated  only  at  the  interface  between  the  two  fluids,  the 

remainder  of  the  flow  remaining  irrotat ional . It  is  convenient 

here  to  integrate  Equation  (11)  across  the  interface  to  yield 

an  expression  for  the  growth  rate  of  the  surface  circulation 

density  a (circulation  per  unit  length  along  the  interface) : 


d a 
dt 


9i°±" 


D 

0 


sin  G 


(12) 


where  o+  is  the  density  to  the  right  of  the  interface  and  o 
that  to  the  left. 


The  total  circulation  of  a given  (ith)  fluid  element 
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. 


dxfdy'is  determined  by 


dr. 

i 


g(p+~  p-) 


dt 


Ay.  n 
1 z 


(13) 


where  Ay^  is  the  length  or  height  of  the  fluid  element  in  the 


vertical  direction  and  nz  is  the'  unit  vector  perpendicular  to 


the  plane  Cf~ motion. 

A convenient  numerical  analysis  of  the  evolution  of  the 
fluid  motion  can  be  obtained  by  dividing  up  the  interface  into 
a number  of  discrete  fluid  elements  and  approximating  the 
circulation  of  each  element  as  being  concentrated  into  a line 


vortex  having  circulation  F.  . The  quantity  Ay.  is  then  to  be 


interpreted  as  the  separation  between  adjacent  vortices.  The 
evaluation  of  the  fluid  motion  then  reduces  to  the  problem 
of  following  the  motion  of  the  individual  discrete  vortices. 

The  velocity  of  the  ith  vortex  is  a summation  over  contributions 
from  all  other  vortices: 


dr . 

l 


N 

V 


dt 


= u . = 

l 


j^i  27r 


F . 

i ,r.  - r . . 

— x (_i j_) 


1 r - r . 
D 


(14) 


This  equation  of  motion  plus  the  relation  determining  the 


circulation  growth  rate  (Equation  13)  in  which  Ay^  is  replaced 


bY  ^yi+i“  yi-i*  yields  a direct  deterministic  procedure  for 
following  the  motion. 


Equations  (13)  and  (14)  have  been  used  to  calculate 
the  evolution  of  a number  of  inviscid,  buoyant  and  shearing 
flows.  These  calculations  are  discussed  in  Section  III. 


10 


a 


H*3  Porous  Medium  or  Viscous  Flows 

When  the  viscosity  forces  dominate  the  inertial  forces 
Equation  (1)  reduces  to 

- -i-  Vp  + nV2  u + g = 0 (15) 

For  the  flow  of  a fluid  between  two  parallel  plates,  the  flow 
is  locally  Poisuelle-like  and  the  viscous  term  is  dominated 
by  the  curvature  of  the  velocity  profile  in  the  direction 
normal  to  the  plates : 

v u * - -T  (16) 

d^ 

where  d is  the  plate  separation  and  u is  the  centerline  velo- 
city. Rewriting  Equation  (15)  we  have 

2 2 

" = ' SWp  Vp  + Hn  5 (17) 


A similar  relation  holds  for  flow  in  a porous  medium.  The 
zero  order  flow  (Vp=Vpo,p=po,n=no)  is  a uniform  velocity  UQ : 

^o  = " 8tTp  Vpo  + ffn  ® {18) 

Taking  the  curl  of  Equation  (17)  we  obtain  for  the  vorticity 


« = I"  [ "v  (??) x v?  + v(^  5J  (19) 

To  first  order  in  the  density  and  viscosity  variations  Vp  may 
be  replaced  by  VpQ  where 

8ri  p 

V?o  = - 3o  + (20) 

d 


X--WL, 


In  other  words 


nasi «;  + £p  x o. 
Vo  0 po  \8n  , 


where  UQ  and  g are  in  the  same  direction. 

For  two  uniform  fluids  separated  by  a sharp  interface, 
the  surface  circulation  density 


and  the  total  circulation  of  a given  (i  ) fluid  element  is 
determined  by 

ri = [(! V1)  u° + (^j  £*]  ^ 1231 


When  the  density  difference  is  not  sufficiently  small, 
such  as  the  case  of  air-sea  interface , to  permit  the  Boussinesq 
approximation  (which  essentially  sets  the  pressure  gradient 
in  the  equation  for  the  oorticity  equal  to  a constant  and 
uniform  value)  a more  complex  procedure  is  required.  We  con- 
sider first  the  somewhat  contradictory  case  of  two  uniform 
density  fluids  separated  by  a sharp  interface  but  having 
negligible  surface  tension.  We  take  the  fluids  to  be  initially 
irrotational . According  to  the  momentum  equation 


- - vi  X vp 


1 


vorticity  will  be  generated  only  where  V — is  non-zero,  i.e., 
at  the  interface.  Also,  according  to  the  momentum  equation, 
the  pressure  gradient  is  partly  due  to  a hydrostatic  head  (as 
in  the  Boussinesq  approximation)  and  partly  due  to  inertial 
effects  in  the  fluid  (not  included  in  the  Boussinesq  approxi- 
mation) : 


Vp 


£>g  - p 


du 

dt 


l 

5 

I 


Thus  the  generation  of  vorticity  has  two  sources:  first,  a 

| 

buoyancy  term  due  to  gravity,  and  second,  an  equivalent 

■jj 

buoyancy  term  due  to  the  fluid  acceleration  (an  "effective" 
gravity)  : 


ri  1 -+  , „ 1 du 

pV  - x g + pV  --  x 
P P 


dt 


= VSnp 


[,-  - «] 
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Explicitly 


This  set  of  equations  (i  = 1 to  N)  is  a set  of  linear  equations 
for  the  unknown  quantities  dr^/dt. 

In  the  limit  of  very  large  density  ratios  ( p 2/Pj_  * °°> 
such  as  for  the  ocean  surface) , the  term  (g  - du/dt)  must 
vanish. 

The  equations  for  the  uorticity  generation  rates 
(diydt)  then  reduce  to  the  density  ratio  independent  result 


(r . -r  .} 
1 J 


In  other  words,  fluid  at  the  surface  slides  freely  over  the 
underlying  fluid  layers  at  a rate  determined  simply  by  gravity 
and  the  local  wave  slope. 


The  momentum 

du  _ _ 1 

eft  p Vp 


equation  has  the  form 
+ g + vV  u + surface  tension  effects. 


In  two  dimensions  for  an  incompressible  fluid,  the  corresponding 
vorticity  equation  is 


az  - 


Vp  + 


surface  tension  effects  + viscous 

diffusion. 


To  estimate  surface  tension  effects,  we  model  the  surface  as 
a layer  of  finite  thickness 


FIGURE  A.  A Finite  Thickness  Interface  between  Two  Fluids 


The  net  surface  tension  force  is  perpendicular  to  the  boundary 
(i.ev,  in  direction  of  density  gradient)  and  is  proportiona1 
to  boundary  curvature  (K) . We  will  model  it  as  a body  force 


acting  within  the  surface  layer  so  that  the  momentum  equation 
takes  the  form 

at  = - F VP  + 5 + vv2S  + vp 

where  T is  the  surface  tension,  K the  local  curvature  (assumed 
to  be  uniform  through  the  layer)  and  Ap  the  density  difference 
(Ap  <<  p)  . To  transform  this  equation  to  vorticity  form,  we 
integrate  along  a line  contour  intersecting  the  boundary  which 
encloses  a fixed  mass. 


U„(s)  U| | (s  + as) 


FIGURE  B.  Variables  Defined  on  a Finite  Thickness  Interface 
Between  Twc  Fluids 

1 | 7+  -*  T f 

— Vp*ds  + v <j)  V u*ds  + <D  KVp*ds 

The  integral  over  the  surface  tension  term  has  contributions 
only  from  the  end  sections  and  reduces  to 
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T 3K 

p Is 


As 


The  viscous  dissipation  term  must  be  treated  carefully. 
Without  the  gravity  and  buoyancy  effects  the  vorticity  equation 
has  the  form 


3§  - 


and  implies  that  any  concentration  of  vorticity  will  spread 

out  by  diffusion.  Thus  the  sharpness  of  the  interface  diffuses 

away  in  time.  In  the  present  model,  we  wish  to  retain  the 

concept  of  a sharp  interface  but  at  the  same  time  to  introduce 

an  effective  viscosity  that  will  ensure  smooth  (or  controllably 

smooth  or  rough)  solutions.  Thus,  for  numerical  purposes  we 

suppress  the  diffusion  in  the  direction  perpendicular  to  the 

2 

interface.  In  other  words,  we  replace  the  term  vV  C by 
2 

g ^ r 

v — *■  where  s is  the  distance  along  the  boundary.  Then, 

as22 

integrating  over  our  mass  elements,  we  obtain  for  the  viscous 
term 


j V2u' 


ds  -*■  v t)  V S*d  (area) 


v 


f 


d2o  As 
ds2 


ds 
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where  is  the  normal  distance  to  the  interface  and  0 is 
the  surface  density  of  circulation  (oAs  = T) . 

The  buoyancy  term,  in  the  Boussinesq  approximation 
(Vp/p  <<  1;  Vp  -*•  -pg)  becomes 


Ap 


g Ay 


where  Ay  is  the  vertical  separation  of  the  ends  of  the  mass 
element  under  consideration.  Thus  the  circulation  of  this 
surface  element  is  given  by 


d r 

dt 


' *p  9 Ay  + v J As  + 1 £r  is 


In  terms  of  the  surface  circulation  density  a,  and  the 
surface  inclination  0 (radians;  tan  0 = dy/ds) 


d_a  _a_  dAs 

dt  As  dt 


2 

Ap  „ da  T dK 

g sin  0 + v — — + - 

P ds2  p ds 


or 


da 

3t 


dAs 


= - a 


3u 


3U|I  Ap  j „ j d2o  T dK 

■3-1.  - _ g Sln  e + - gj 


Since  - As  ^ — and  since  on  the  surface’ 


d 

dt 


at 


+ Un 


3s 


For  finite  a there  is  a discontinuity  across  the  surface  in  the 
parallel  velocity.  The  value  for  up  in  these  equations  is  the  mean 

[uj|  = i(Ujj  ( + ) + u..  (-))]  . 
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we  may  write 


3a 

St 


3a  3ul  Ap  „ 32a  T 3K 

_U||  3l'flSr'"p9sin0+v-T+-Ti. 


3s 


(ou„)  - A£ 


3 o T 3K 

Si"  vuu||  i - 7 9 sin  0 + v - 7+_T_  . 

o S 


Previously  we  have  followed  specific  mass  elements  in 
their  motion  (Lagrangian  formulation) . It  is  also  convenient 
to  divide  the  surface  infcd  fixed  length  arc  segments  (mixed 
Eulerian-Lagrangian) . Let  sq(s)  be  the  original  arc-length 
of  the  mass  Element  which  is  now  at  the  arc-length  s measured 
fror  a stagnation  point  (i.e.,  a point  where  uB  = 0 for  all 

5 

time).  We  want  to  develop  equations  for  two  quantities: 

1)  the  original  arc-length  as  a function  of  the  present  arc- 
length  s (s,t),  and  2)  the  surface  coordinates  r(s,t).  For 
the  function  sQ(s,t),  we  note  that  if  we  write  s t s(s  ,t), 
then  for  fixed  s 


ds  = 0 


Thus* 


(3s/3t) 

so 

(3s/3s  ) 


(“I,  «*o>  - un  < 


* 


Subscripts  identify  the  variable  being  held  fixed  during 
the  partial  differentiation. 


where  Ujj  (sQ)  is  the  parallel  velocity  at  the  point  r(sQ) 
(although  Uy (0)  is  0 according  to  our  stagnation  point 
reference,  we  keep  the  general  form  here).  Also  we  may 
write  for  r = r(s,t) 


Thus,  the  three  equations  (taking  u.,  (0)  = 0) 


3r 

at  j 


= u - u 


3r 


II  9a  L 


3o_ 

at, 


- (ou„  ) - g sin  0 + v — a-  + - 9K 


3s  w“|| 


3s 


2 p 3s 


u(r,t)  = 


= f t (s,t)  x 

J I r - r'| 


ds 


permit  following  the  evolution  of  the  interface  in  time.  In 
addition,  the  relation 


may  be  integrated  to  indicate  the  degree  of  flow  movement  along 
the  boundary. 


HI  SCALING  AND  LABORATORY  SIMULATION 

Vortex  consideration  furnishes  a powerful  ally  in  attacking 
many  of  the  complex  problems  of  non-linear  rotational  flows, 

We  shall,  in  this  context,  establish  a working  and  efficient 
numerical  basis  for  such  approach  by  emphasizing  the  manner 
in  which  flow  motion  is  generated  by  the  vorticity  and  how 
the  subsequent  evolution  develops. 

In  this  section  we  study  numerically  models  of  the 
flows  discussed  in  the  introduction.  In  Section  III.i,  the 
evolution  of  an  injection  cylinder  is  shown  to  result  in  a pair 
of  line  vortices,  in  Section  III. 2,  the  buoyant  rise  of  a 
cylinder  of  heated  gas  is  shown  to  result  also  in  the  develop- 
ment of  a line  vortex  pair.  In  Section  III. 3,  the  Kelvin- 
Helmholtz  vortices  of  waves  generated  in  a shear  layer  are 
calculated  for  three  different  Richardson  numbers.  In  III. 4, 
we  repeat  the  original  calculation  of  that  of  Rosenhead  (1931) 
and  that  of  Kronauer  and  Abernathy  (1962)  showing  the  develop- 
ment of  a Karman  vortex  street,  in  Section  III. 5,  the  well 
known  development  of  tip  vortices  is  simulated.  In  III. 6, 
we  study  a rather  unique  flow  in  which  a cylinder  of  intermediate 
density  fluid  is  placed  at  a stable  interface  between  two  fluids. 
The  collapse  of  this  cylinder  as  it  seeks  its  own  level  results 
in  a splitting  of  the  original  cylinder  into  two  laterally- 
moving  element,  each  of  which  consists  essentially  of  a line 
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vortex  pair.  In  Section  III. 7,  a study  of  the  same  cylinder 
as  that  in  III. 7 is  shown,  but  it  is  placed  at  some  distance 
below  the  thermocline.  In  Section  III.  8,  finally,  the  develop- 
ment of  Taylor-Saf fman  instabilities  at  an  interface  between 
two  viscous  fluids  in  creeping  flow  is  demonstrated. 

Rise  of  an  Injection  Cylinder 

The  formation  of  a ring  or  a pair  of  line  vortices 
following  a sudden  introduction  of  a blob  of  fluid  into  a 
quiescent  medium  of  the  same  density  or  an  impulse  given  to 
the  surrounding  fluid  is  a commonly-observed  phenomenon. 
Examples  of  these  are:  a smoke  ring,  a pulsating  jet,  or  a 

passage  of  a high-speed  streamlined  vehicle.  In  this  last 
example,  a column  of  ambient  air  will,  in  addition  to  the 
axial  motion,  be  pushed  into  vertical  ascent.  The  ensuing 
motion  will  be  dominated  by  the  vorticity  generated  as  a result 
of  the  impulsive  shearing  motion.  To  simplify  the  problem, 
one  could  conceive  a cylinder  of  air  impulsively-injected 
upward  so  that  the  cylinder  has  a uniform  velocity.  If  one 
fixes  the  coordinate  on  the  cylinder,  the  external  flow  is 
simply  a potential  flow  past  a circular  cylinder  of  which  a 
solution  is  given  as  a result  of  a vortex  doublet  (Batchelor, 
1967,  p.  535).  Since  in  inviscid  two-dimensional  and  non- 
buoyant  flow,  the  vorticity  is  constant  throughout  the  motion, 
it  is  convenient  to  attach  the  vorticity  to  the  shear  interface 
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in  order  to  trace  the  boundary  of  the  cylinder.  A distribution 
of  the  vorticity  along  the  boundary  can  be  found  to  yield  a 
uniform  velocity  in  the  cylinder.  We  shall  return  to  this 
point  after  we  have  non-dimensionalized  the  governing  equations. 

Figure  1.1a  shows  the  initial  distribution  of  the 
vortices,  with  the  vortex  axes  being  parallel  to  the  generatrices 
of  the  cylinder.  By  symmetry,  the  vorticity  is  of  the  opposite 
signs  on  each  half  of  the  circle.  Since  there  is  no  density 
difference  involved  in  this  case , the  vorticity  is  constant 
throughout  the  motion.  The  motion  is  determined  by 

either  Equation  (4)  or  by  the  set  containing  Equations  (6  ) 
and  ( 7 ) . These  equations  are  reduced  to  dimensionless  form 
by  introducing  the  following  characteristic  dimensions: 

distance  R:  initial  Cylinder  radius 

time  T : 

1 o 

circulation  rQ:  total  circulation  assumed  to  be 

distributed  on  the  half  circle. 
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In  terms  of  the  dimensionless  distance  £ = n = the 

K K 

dimensionless  time  t = ^ and  dimensionless  circulation 
r 

Y = — / then  Equation  ( 4 ) becomes 


dr 


. - I 


Yj  ( n j ~ni ) / [ (5j“5i)  2 + (ilj-r^)2] 


(24) 


dn. 


dx 


z 


Yj(^j-?i)/[(^j-?i)2  + (nj-r^)2]  (25) 


for  the  Green's  function  formalism.  In  terms  of  the  dimension- 


2nT 


less  stream  function  Y = ^_L  and  the  vorticity  C = — , 

o R T 

the  Equations  ( 6 ) and  ( 7 ) become  ° 


dq 

~ } 
94' 

(26) 

dx 

9n  i 

i 

dni 

9y| 

(27) 

dx 

h] 

i 

2~ 

V V = 

l 

(28) 

for  the  stream  function  formalism. 
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The  cylinder  boundary  at  time  zero  was  divided  uniformly 
into  N points  (61  pts.  over  the  half  circle).  Each  point  is 
assumed  to  have  a constant  vorticity  according  to  * j( ^i+l”71!-!^  ' 

i = 2,  N - 1,  and  y^  and  yN  are  ec3ual  to  zero.  From  this  dis- 
tribution of  discrete  vortices,  the  velocity  of  each  vortex  can 
be  calculated  using  Equations  (24)  and  (25) . To  verify  the 
uniformity  of  the  velocity  inside  the  cylinder,  Figure  1.1b 
shows  the  initial  velocity  distribution.  The  positions  of 
vortices  are  then  advanced  through  the  integration  of  the 
obtained  velocities.  If  stream  function  is  desirable,  it  can 
be  calculated  using  Equation  (28)  , and  from  ¥ the  velocity  can 
be  obtained  for  the  next  time  step. 

The  vortices  tend  to  be  coagulated  to  the  bottom  center 

of  the  cylinder  initially,  and  the  purely  upward  translation 

does  not  take  place  as  it  occurs  in  the  thermal  which  will  be 

discussed  in  the  next  section.  At  t = .2,  the  cylinder  evolves 

into  the  well-known  kidney  bean  shape.  This  is  shown  in 

Figure  1.2.  In  Figure  1.3,  the  center  of  vorticity  is  well 

developed  at  t * .4  and  in  Figure  1.4  the  vortices  have  rotated 

around  the  center  of  the  vorticity  several  times  at  t = 1. 

Batchelor  (1969)  showed  that  the  centroid  of  the  vorfciilty 

should  be  a constant  of  motion.  In  this  study,  it  appears  that 

— 2R  — 

the  centroid  of  half  of  the  circle  x = — , y = 0 indeed  remains 
to  be  constant. 

We  have  utilized  vortices  of  various  finite  core  radii 
to  obtain  the  velocity  at  each  vortex  from  the  Green's  function 
formalism  in  Equations  (24)  and  (25)  . When  the  core  radius  is 
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increased,  the  net  effect  on  the  flow  is  found  small  on  the 
large  scale  but  the  flow  in  the  small  scale  does  become 
smoother . 


III. 2 Rise  of  a Buoyant  Cylinder  (Thermal) 
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FIGURE  1.1 


Initial  Vortex  Position  of  an  Injection 
Cylinder 
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FIGURE  1.2  Injection  Cylinder  at  T 
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floating  upwards  under  the  influence  of  gravity.  The  fluid 
within  the  cylinder  is  assumed  initially  to  have  a uniform 
density  (p^)  and  to  be  imbedded  in  a slightly  more  dense 
fluid  (also  of  uniform  density  p2> . The  fluids  are  assumed 
to  be  miscible  (no  surface  tension) . In  this  case,  Equations 
(13)  and  (14)  may  be  reduced  to  dimensionless  form  by  intro- 
ducing the  following  characteristic  dimensions: 

distance:  R = the  inital  cylinder  radius 

time:  T = [g  (Ap/p) /2ttR] 

circulation:  Fq  = [2'nq  ( Ap/p)  R3  ] 1^2 

In  terms  of  the  dimensionless  distances,  £ = x/R,  n = y/R, 
the  dimensionless  time  T = t/T,  and  dimensionless  circulation 
y = r/rQ,  the  equations  of  motion  (13)  and  (14)  become 


(29) 


)/[U.-ZL)2  + (Hj-n^2]  (30) 


^/(Uj-q)2  + (nj-r^)2]  (31) 


where 


Equations  (29)  to  (31)  have  been  used  to  calculate  the 
time  dependent  motion  in  two  dimensions  following  the  release 
of  an  initially  uniform  circular  cylinder  of  light  weight  fluid 
in  a homogeneous  heavier  fluid.  Since  the  density  gradients 
in  this  example  are  limited  to  the  (deforming)  surface  of  the 
cylinder,  the  motion  may  be  followed  by  following  the  history 
of  the  vortex  sheet  which  comprises  the  cylinder  boundary  (see 
Fig.  2.1).  The  results  of  the  calculation  are  shown  in 
Figures  2.2  to  2.5.  The  cylinder  boundary  at  time  zero  was 
divided  uniformly  into  N points  (61  points  on  the  half  circle) . 
The  velocity  of  each  point  was  calculated  at  successive  time 
increments  according  to  Equations  (30)  and  (31) . The  circula- 
tion of  each  point  was  calculated  from  Equation  (29)  as  a 
function  of  time,  the  initial  values  being  taken  equal  to 
zero  (i.e.,  no  initial  motion). 

The  initial  motion  of  the  cylinder  appears  to  be  simply 
an  upward  displacement  without  sensible  distortion.  By  the 
time  the  net  displacement  is  of  the  order  of  1/2  the  initial 
cylinder  radius,  the  beginning  of  vortex  development  is  evi- 
dent (Figure  2.3).  The  vortex  appears  well  developed  by  the 
time  the  buoyant  region  has  risen  about  one  diameter  (Figure 
2.4,.  By  this  ‘ ime  most  of  the  vorticity  is  concentrated  in 
the  vortex  region.  The  rate  of  change  of  the  total  circula- 
tion of  this  region  is  obtained  by  summing  Equation  (29)  over 
the  entire  vortex  sheet 


33 


^(?Yi)  = 5n 

or  (32) 


where  6n  is  the  thickness  of  the  cap  on  the  axis  of  symmetry. 
Thus,  the  value  of  the  vortex  circulation  grows  during  the 
vortex  development,  but  saturates  when  the  vortex  has  fully 
developed.  The  subsequent  motion  (after  vortex  formation) 
of  the  buoyant  and  entrained  material  has  been  discussed  by 
others  [particularly  by  J.  S.  Turner  (1959)  and  also  by  T. 

Fohl  (1967)]. 

Although  the  present  calculation  was  carried  out  for 
a cylindrical  configuration  essentially  similar  results  are 
anticipated  for  spherical  buoyant  bubbles. 

The  entrainment  process  involved  in  this  simple  inviscid 
model  thus  appears  to  be  a simple  enfolding  of  the  ambient 
fluid.  Turbulence  effects  may  alter  the  processes  somewhat f 
particularly  in  determining  the  detailed  structure  and  degree 
of  mixing  within  the  vortex  core.  Careful  experiments  are 
valuable  for  developing  models  that  include  these  effects. 

Both  experimental  and  numerical  studies  should  be  performed  to 
determine  the  effects  of  large  initial  density  differences, 
finite  initial  density  gradients  in  the  bubble,  atmospheric 
stratification  and  wind  shear,  finite  initial  turbulence, 
and  finite  initial  translational  velocity. 
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At  times  later  than  2T,  when  the  cylinder  (see  Fig.  2.5 
has  risen  more  than  one  ciameter  the  set  of  point  vortices 
form  an  irregular  distribution  within  a finite  cloud.  The 
original  vortex  sheet  is  now  so  convoluted  as  to  be  impossible 
to  follow.  Although  the  numerical  model  cannot  be  a good 
model  of  the  small  scale  structure  at  such  times  it  is 
interesting  to  note  that  the  large  scale  motion  agrees 
reasonably  well  with  theoretical  expectations  (at  least  to 
values  of  t < 3T) . This  may  be  seen  as  follows. 

Turner  (1957)  has  shown  that  the  circulation  of  each 
vortex  approaches  a constant  value  after  vortex  formation. 

This  may  be  seen  from  Kelvin's  theorem  which  states  that 
around  any  closed  .circuit  C 

dr  f 1 n ^ 

at  = f p '7p-ds 

J c 

After  vortex  formation,  the  density  along  a path  threading 
the  center  of  the  vortex  is  essentially  constant  and  equal 
to  the  ambient  value  and  dT/dt  -*•  0.  When  F is  constant,  the 
rise  velocity  varies  inversely  as  the  separation  of  the 
vortex  pair 

V ~ 1/R 

The  upward  momentum  increases  at  a constant  rate 


where  Fg  is  the  (constant)  buoyant  force. 

2 

Since  M is  proportional  to  R in  two  dimensions  and  RV  is 
constant,  the  separation  R increases  linearly  with  time 

R ~ t . 

Since  the  rise  velocity  of  the  vortex  pair  varies  as  F/R, 
the  net  rise  distance  y increases  logarithmically  (in  two 
dimensions)  with  time 

y ~ 9.n  t 

In  Fig.  2.6,  we  show  that  the  time  dependence  of  the  width 
and  height  of  the  rising  vortex  pair  agree  reasonably  well 
with  the  expected  dependence. 

In  three  dimensions  the  expansion  rate  will  have  a 
different  time  dependence.  Since  here  the  mass  varies  as  R 
the  momentum  equation  reduces  to 


after  torus  formation  (when  RV  ~ constant) . Here  R ~ t1^2 
-1/2 

and  dz/dt  ~ t . Thus  the  radius  of  the  torus  increases 
linearly  with  height  (R  ~ z)  . 

When  the  fast  Fourier  transform  (FFT)  is  applied  to 
solve  the  stream  function,  a vortex  system  of  much  larger 
numbers  of  particles  can  be  employed  costing  essentially  the 
same  amount  of  computation  time.  For  example,  in  this  case, 


the  Green  s function  approach  using  61  points  took  .246 
seconds  per  time  step.  The  stream  function  approach  using 
591  points  took  .426  seconds  per  step,  while  the  same  approach 
using  4.1  points  costed  .31  seconds  per  step.  One  drawback 
is,  however,  that  the  finite  mode  Fourier  transform  does 
produce  aliasing  errors.  Figures  2.7  through  2.11  show 
the  result  obtained  from  the  stream  function  using  200  points 
through  the  same  period  as  that  in  Figures  2.1  to  2.5. 

The  velocity  vector  plot  seems  to  be  smooth  but  the  vortex 
position  plot  shows  there  are  small  fluctuations  developing. . 
If  they  are  not  eliminated,  these  fluctuations  will  be 
amplified  into  large  amplitude  errors.  We  applied  a smoothing 
function  to  high  wavenumber  portions  of  the  Fourier  components 
to  eliminate  this  noise  whic.i  is  well  known  as  the  Gibb's 
phenomenon.  The  result  of  the  damping  is  a smoother  rolled 
up  configuration. 
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FIGURE  2.4  Vortex  Pair  Developed  from  a Buoyant 
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FIGURE  2.8b.  Velocity  Vector  Plot  of  a Buoyant 
Cylinder  at  x = .72 


ISCRETE  VORTEX 


'^v'wciV  vVto%  'l*t  ♦ ♦A4****** 

AT  TH8  T I MB ■ I . • t 0 0 

► Auf*8CAXB  L ■ H 0 T H imiONTI  T H ■ MAXIMUM  FLOW  BP«b6  lITUMrtB  W ^ 

^ ^ ii^mt  .{ftafeiMfe  ok  im^cia^ula4io^li^tmA  •Al^o  ^ ^ ^ ^ 


\ \ \ t t f / / / 


^ K \ ^ \ { f t 

* * * ; ; It  I t 1 1 4 ; ; ; ? 

* ♦ * * A f t '■  '■  * * 1 

: i t tuN.,:' 


^ * 

^ * 

T * A 4,  1 


0 


1 1 1 • 3 Finite  Amplitude  Kelvin-Helmholtz  Waves 

As  early  as  1868,  Helmholtz  found  that  the  vortex  sheet 
of  infinite  horizontal  dimension  formed  by  shear  is  unstable 
to  any  kind  of  disturbances.  Recent  studies  of  the  mountain 
lee  waves  and  clear  air  turbulence  have  renewed  interest  in 
studying  the  evolution  of  Kelvin-Helmholtz  waves  in  strati- 
fied mediums.  Radar  backsca ttering  studies  have  revealed 
the  formation  and  growth  of  the  cat's  eyes  in  the  atmosphere 
(Richter,  1969)  and  Kelvin-Helmholtz  billows  are  also  found 
by  Wood  (1969)  on  the  seasonal  themocline  in  the  Mediterranean 
Sea  off  the  coast  of  Malta.  A numerical  study  of  Kelvin- 
Helmholtz  waves  in  a viscous  luid  which  solves  the  Navier- 
Stokes  equations  was  given  by  Patnaik  (1973) . 

In  this  section,  we  present  the  results  of  a study 

of  this  subject  in  an  inviscid  fluid  at  three  effective 

Richardson  numbers,  -1,  0,  100,  which  are  defined  as  R = - 9APr3 
n /r.4  \ i . „ 2 


The  equations  of  motion  are  based  on  Equations  (24)  and 
(25)  for  the  case  R..  = 0,  and  Equations  (29),  (30)  and  (31) 
for  the  case  when  R^  f 0 . The  infinite  vortex  sheet  is  now 
subjected  to  periodic  disturbance  in  vertical  displacement. 

The  sum  of  the  velocity  contribution  from  all  the  vortices 
located  at  multiple  wavelengths  R,  2R,  3R,  ...  apart  is  evaluated 
by  the  infinite  series  result  given  by  Lamb  (1932) . The  vortex 
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sheet  is  divided  uniformly  into  41  points  and  each  contains 
initially  an  equal  amount  of  vorticity  generated  by  the  shear 
across  it. 

(i)  R^  = 0 case  (pure  shear  case) 


The  vortex  sheet  is  given  an  initial  disturbance 
shown  in  Figure  3.1,  a sine  wave  of  wave  length  R = 2,  and 
amplitude  §•  = .1. 

K 


When  the  vortex  sheet  is  rolling  up,  it  reaches  the 
breaking  height  at  i = .64;  that  is,  t = .641^ — I where 


Fq  is  the  total  circulation  around  one  R length  of  the  vortex 
sheet.  This  is  shown  in  Figure  3.2.  At  t = 1.04,  as  shown  in 
Figure  3.3,  one  complete  turn  is  made,  and  Figure  3. 4 shows 
the  final  rolled  up  configuration  at  ' = 2. 


(ii,  R^  = - 1 (unstable)  case 

This  corresponds  to  the  case  when  the  fluid  is  unstably 
stratified;  therefore,  large  amplification  of  the  disturbance 
and  rapid  roll  up  should  occur.  The  same  initial  disturbance 
is  given  to  the  vortex  sheet  as  in  the  R^  = 0 case  (Fig.  3.5)  , 
that  the  fluid  is  stratified  with  the  lower  density  fluid 
underlying  the  heavier  fluid.  Since  R^  < 0 , the  Brunt -Vaisal; 
frequency  becomes  imaginary  so  that  no  oscillatory  motion  can 
exist.  Figure  3.6  shows  that  the  vortex  sheet  reaches  the 
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except 


breaking  point  at  x - -57,  which  is  earlier  than  t = .64 
for  the  = 0 case.  Figure  3.7  indicates  that  at  t = .81, 
one  over  turn  has  been  completed.  Figure  3 . 8 shows  that 
the  wave  has  grown  into  a much  larger  roll  at  t = 2 than 
that  shown  in  Figure  3.4. 

(iii)  = 100  (stable)  case 

This  corresponds  to  the  case  when  the  atmosphere  is 

very  stably  stratified,  but  with  strong  wind  shear.  It  is  an 

ideal  situation  for  the  generation  of  the  internal  gravity 

waves.  One  will  not  find  any  roll  up,  only  oscillatory  motion 

at  fixed  Brunt-Vaisala  N =-yj~  The  frequency  N is  easily 

determined  from  the  definition  of  the  Richardson  number, 
r-  ro 

that  is,  l'  = R — j ; hence,  it  is  expected  to  take  t 4 
' R^ 

to  complete  one  cycle  of  oscillation. 

Again,  the  same  initial  disturbance  as  shown  in  Figure  3.1 
is  imposed  upon  the  vortex  sheet  which  now  possesses  only  ten 
percent  of  the  circulation  of  either  of  the  previous  two  cases. 

At  t = .99,  that  is  at  one  quarter  of  the  period  of  the  oscilla- 
tion, the  vortex  sheet,  which  is  shown  in  Figure  3.9,  becomes 
simply  a straight  line,  and  at  x = 2.01  (Figure  3.10),  half  the 
cycle  of  the  oscillation  is  completed. 

So  far  the  calculation  includes  only  one  disturbance  wave- 
length. For  practical  applications,  it  may  be  of  interest  to 
calculate  interactions  of  disturbances  with  different  wavelengths  in 
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FIGURE  3.2 


Kelv in-Helmholtz  Wave  in  Pure  Shear 
Case  at  f = .64 
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FIGURE  3.9  Kelvin-Helmholtz  Wave  in  a Stably 
Stratified  Medium  at  x = .99 
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order  to  find  the  fastest  growing  mode.  The  mountain  lee 
wave  problem  can  also  be  treated  conveniently  through  the 
present  approach. 


HI. 4 Development  of  the  Karman  Vortex  Street 

It  is  well  known  that  the  vortex  wake  of  a circular 
cylinder  becomes  unsteady  once  the  Reynolds  number  based  upon 
the  cylinder  diameter  exceeds  40  and  vortex  street  will 
develop  subsequently.  Karman  (1911)  pointed  out  that  whan 
the  Reynolds  number  increases  (hence  the  viscosity  effects 
decrease)  the  vorticity  contained  in  each  of  the  wake  vortices 
will  not  be  dissipated  rapidly  enough  to  prevent  the  vortex 
interaction  among  themselves.  That  is,  for  high  Reynolds 
number  flows,  the  study  of  the  viscous  wake  can  be  treated 
as  in  an  inviscid  fluid  with  the  motion  dominated  by  the 
vortex  interactions  between  the  two  vortex  sheets.  These 
sheets  of  opposite  signs  shed  from  the  moving  body  in  a viscous 
fluid  can  be  defined  as  inviscid  vortex  sheets  located  at  the 
velocity  inflexion  points.  Since  that  is  where  it  contains 
the  most  of  the  vorticity,  the  vortex  sheets,  once  under  the 
disturbances  of  different  wavelengths,  will  evolve  into  different 
vortex  street  configurations.  Among  them,  the  most  stable  one 
- and  therefore  the  most  often  observed  one  - the  Karman  vortex 
street  should  have  a ratio  of  the  separate  distance  h to  the 

u. 

disturbance  wavelength  A , j to  be  equal  to  .281. 
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The  fact  that  the  two  vortex  sheets  are  unstable  to 
any  disturbances  was  also  predicted  by  the  Orr-Sommerfeld 
stability  theory  which  states  that  any  velocity  profile  having 
an  inflexion  point  is  unstable,  so  that  an  unsteady  analysis 
based  upon  the  central  theme  of  the  vortex  interactions  should 
be  carried  out.  The  often  observed  unsteady  separation  bubble 
at  a concave  corner  is  a good  example  in  this  context. 

Figure  4 .1  shows  the  initial  configuration  of  the  two 
vortex  sheets  displaced  vertically  by  a sine  wave  of  wave 
length  h = .28  A;  h is  the  separation  distance  between  the 
two  sheets;  the  wave  amplitude  is  ^ = .1.  The  same  summation 
given  by  Lamb  (1932)  as  mentioned  in  Section  III. 3 is  utilized 
in  this  section  also.  Figure  4.2  shows  that  at  x = .62, 

Kelvin -Helmholtz  waves  reach  the  breaking  height.  The  dots 
and  triangles  are  used  to  distinguish  the  signs  of  the  vortex. 
Figure  4.3  shows  that  at  x = 1.2  most  of  the  vorticity  is 
coagulated  into  the  al ternatingly-space  horseshoe  vortices. 

Figure  4.4  shows  the  vortex  street  formed  at  x = 2.10;  the 
spacing  between  the  vortex  clouds  becomes  clearer  and  we  notice 
that  the  vortices  of  the  opposite  sign  are  mingled.  There  is 

r 

also  a pure  translation  which  is  determined  by  =2.  toward  the  left 
Abernathy  and  Kronauer  (1962)  found  noisy  results  when 
using  vortices  of  zero  core  radius.  We  have  obtained  smoother 
roll  up  by  using  vortices  of  finite  core  radius.  The  core 
radius  is  not  a critical  parameter  in  the  determination  of  the 
solution.  We  have  conducted  tests  using  different  core  radii 
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and  found  that  the  overall  motion  of  the  vortex  system  is 
independent  of  the  core  radius  but  the  small  scale  motion 
loses  its  random  character  if  core  radius  is  large.  Physical 
argument  is  required,  however,  to  determine  the  solution  at 
small  scale.  It  is  noted  that  unless  a three  dimensional 
viscous  VIC  method  is  used,  the  solution  obtained  at  small 
scale  is  not  meaningful.  By  choosing  a core  radius  too  large, 
the  solution  loses  its  accuracy.  But  by  choosing  too  small  a 
core  radius,  the  solution  will  become  noisy.  Although  the 
two  dimensional  VIC  fails  to  predict  the  small  scale  motion, 
one  can  still  get  a fascinating  glimpse  into  the  randomization 
from  the  nearly  random  small  scale  structure  within  an  organized 
large  scale  structure  as  a result  of  vortex  interaction  (as 
pointed  out  by  Liepmann,  1961)  and  also  a demonstration  of  the 
statistical  characteristics  of  the  turbulence  proper.  That  is, 
the  energy  is  transferred  from  large  scale  to  small  scale. 

Figure  4.5  shows  the  initial  configuration  for  the  case 
h 

when  X = Figs*  4.6  and  4.7  show  the  solution  at  t=.92,  1.5 

Figure  4.  8 shows  the  result  at  t = 2.14.  Notice  that  the  vortices 
have  not  yet  developed  into  the  stable  configuration  as  shown 
by  Figure  4.4.  Further  coagulation  will  take  place  and  will 
eventually  lead  to  the  Karman  street. 

The  fact  that  the  VIC  method  predicts  the  lateral 
broader ing  of  the  vortex  street  in  the  absence  of  viscosity 
and  turbulence  is  a promising  feature  for  the  turbulence 
modelling  by  the  VIC  method. 
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FIGURE  4.4  Distribution  of  Vortices  in  a Karmen  Vortex 
Street 
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FIGURE  4.5 


Initial  Disturbance  j-  = .12  on  a Pair 
of  Vortex  Street  A 
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FIGURE  4.8  Development  of  Vortex  Street  for 


III. 5 Aircraft  Trailing  Vortex  Street 


By  Prandtl's  lifting  line  theory,  the  trailing  vortex 
shed  from  the  low  aspect  ratio  wings  are  straight  and 
parallel  to  the  direction  of  flight  and  the  flow  in 
the  neighborhood  of  any  one  section  of  the  wing  is  approxi- 
mately two  dimensional  and  independent  of  the  neighboring 
sections.  This  seems  to  be  valid  for  most  of  the  present 
day  transport  aircraft  (excluding  SST's);  however,  for  the 
low  aspect  ratio  wings,  the  framework  laid  out  in  tnis 
study  is  still  applicable  except  that,  a three-dimensional 
formulation  must  be  utilized.  Until  practical  application 
warrants  the  complication,  we  shall  assume  the  hypotheses 
in  Prandtl's  theory  applies. 

By  wing  theory,  the  lift  or  the  wing  loading  is 
linearly  proportional  to  the  circulation  about  the  wing 
cross  section  and  it  is  well  known  that  the  wing  load  [so 
is  the  circulation  S(x) ] can  be  approximated  by  the 
elliptic  curve 


where  SQ  is  the  maximum  circulation  at  x = 0 and 
R is  the  wing  span.  By  Stokes'  law  of  conserva- 
tion of  the  circulation,  the  circulation  at 


x+dx  is  decreased  by  the  amount  As  = S(x+dx)  - S(x)  so  that 
this  amount  of  circulation  must  be  shed  from  the  wing 
section  between  x and  that  at  x+dx  in  the  form  of  cylindrical 
tip  vortex.  Section  AA'  in  Figure  5.1  will  have  a vortex 
sheet  with  variable  strength  resulting  from  the  vortex 
shed  from  the  wing, and  the  net  effect  of  this  vortex  sheet 
is  to  generate  a discontinuity  in  the  horizontal  velocity  u 
along  the  wing  surface. 

Such  vortex  sheet  is  the  trailing  vortex  commonly 
referred  to  as  the  aircraft  wake  vortex;  this  should  not 
be  confused  with  the  trailing  vortex  shed  from  a two- 
dimensional  wing  section  due  to  Joukowski's  analytic 
condition.  Section  BB'  in  Figure  1 shows  this  trailing 
vortex  sheet.  It  is  located  at  z — °°  and  its  strength 
has  the  same  distribution  as  the  wing  circulation  S(x) 
but  of  the  opposite  sign  to  balance  the  infinite  tan- 
gential velocity  at  the  trailing  edge  of  the  wing  flaps. 

From  this  argument  it  is  clear  that  the  circulation  of 
the  trailing  vortex  at  Section  AA'  is  equal  to  the  rate 
of  change  of  S(x) , that  is  - — — — - with  the  proper  sign. 
Direct  differentiation  oi  Eq.  (13)  with  respect  to  the  x 
will  lead  to  infinity  at  the  wing  tip,  hence  the  following 
relation  is  applied  to  obtain  finii-  > . . iation  strength 


I’(x.  ) 


S 


(x. 


i + 1 


( 34) 
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The  vortex  sheet  is  divided  into  NTR  strips  along 
the  z direction,  each  segment  of  which  contains  a circu- 
lation r\(x)  given  by  Eq.  (34). 

The  wind  profile  near  the  ground  is  known  to  ex- 
hibit a logarithmic  dependence  upon  the  elevation  y 
(Blackadar  and  Tennekes,  1968) 


u 


U = 


Jin 


+ constant 


( 35) 


where  ut  is  the  friction  velocity  and  is  usually  given  by 


the  relation  u s — 


U 


, the  k is  the 


t 30  at  height  of  1 km 
Karman  constant  and  is  equal  to. 42  for  most  applications 

(Hinze,  1959) . The  yQ  is  the  roughness  parameter  and 

for  typical  atmospheric  conditions  is  about  .01  m. 

From  Eq.  (35)  the  vertical  wind  shear  can  be  obtained 

by  taking  the  derivative  with  respect  to  y which  yields 


r (y)  = 


u 


<y 


(36) 


This  circulation  is  assigned  to  each  mesh  point  as  a dis 


. 


? 


Crete  vortex.  Notice  the  vorticity  is  infinite  at  y = 0* 


in  order  to  avoid  this  we  applied  the  following  relation 
instead. 


r (o) 


and 


for  j = 1 


r (y-j) 


Ay 

Ay 


\ 


j = 2, 17  (37) 


A study  of  the  trailing  vortex  shed  from  a Boeing 
747  aircraft  at  a height  of  61  meters  (200  ft.)  above  the 
runway  using  a 32x32  grid  was  carried  out.  The  trailing 
vortex  was  represented  by  25  discrete  vortices  over  half 
of  the  wing  span,  each  assinged  a circulation  value 
according  to  Eq.  (34).  The  wind  shear  vorticity  is 
distributed  over  the  flow  domain  on  a 17x32  mesh,  and 
the  images  are  obtained  by  the  symmetry  condition  in  the 
vertical  direction.  The  buoyant  engine  exhaust  is 
also  represented  by  25  vortices;  the  temperature  difference 
with  respect  to  the  ambient  air  is  assumed  to  be  10°K. 
Figure  5.2a  shows  a Boeing  747  trailing  vortex 
and  its  buoyant  exhaust  in  a 17-32  grid. 
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On  each  grid  point  there  is  a wind  shear  vortex  with 
strength  determined  by  Eq.  (37).  The  four  downward 
arrows  indicate  the  reference  points  of  the  initial 
yeometry ; all  the  dimensions  are  in  kms  units.  Figure  5.2b 
shows  the  velocity  vector  plot  including  the  trailing 
vortices,  wind  shear  vortices  and  their  images;  the  maxi- 
mum flow  speed  is  represented  by  the  length  indicated  on 
the  upper  left  corner.  From  this  plot,  it  is  clear  that 
the  ground  does  possess  some  translation  and  the  wind  is 
blowing  from  the  right  to  the  left.  There  is  a vertical 
downwash  induced  by  the  lift  on  the  wing  and  the  wind 
profile  is  significantly  altered  by  the  presence  of  the 
trailing  vortices  - notice  the  flow  is  opposite  to  the 
wind  direction  under  the  upwind  tip  vortex.  Figure  5.3 
shows  the  rolling  up  of  the  vortex  sheet  after  2.08 
seconds,  no  skewness  is  observed  at  this  time,  and  the 
exhaust  plumes  are  elongated  along  the  trailing  vortices. 
Figure  5.4  shows  the  overall  picture  of  the  vortex  system 
at  - 4*08  seconds,  the  wind  shear  vortices  near  the 
ground  where  the  vorticity  is  maximum  are  swept  up  and 
mutual  induction  between  those  wind  shear  vortices  and  the 
tip  vortices  may  be  expected  to  emerge.  Figure  5.5a 
shows  the  skewed  configuration  at  t = 10  seconds  and 
the  exhaust  plume  is  completely  wrapped  into  the  tip 
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vortices . The  trailing  vortices  are  transported  nearly 
40  meters  to  the  left  from  the  original  position,  and 
the  position  of  the  wind  shear  vortices  delineate  clearly 
the  wind  profile.  Figure  5. 5b  shows  the  velocity  vector  plot 
at  this  time.  Notice  that  if  the  vortex  system  is  swept 
out  of  the  boundary  of  the  flow  domain,  the  periodic  con- 
dition implied  by  the  Fourier  transform  will  require  the 
vortices  to  be  replenished  into  the  domain  but  at  one 
periodic  length  apart  from  the  original  position.  One 
can  always  choose  a domain  large  enough  to  avoid  the  in- 
fluence of  the  periodic  images. 
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FIGURE  5.5b.  Aircraft  Lrailinq  vortex  configuration 


Collapsing  Wake  on  an  Ocean  Thermocline 


We  now  turn  to  a phenomenon  which  is  commonly  treated 
as  a hydraulic  jump  problem  rather  than  from  the  point  of  view 
of  the  vortex  interactions.  It  corresponds  to  the  intrusion 
of  a heavier  fluid  (a  front  or  'nose')  into  a fluid  of  lighter 
density.  Examples  of  this  flow  are  found  in  the  atmosphere, 
in  a weather  front  (say,  a sea  breeze) , in  front  of  a gravity 
current  which  is  usually  termed  "Sudanese  habcob",  at  river- 
sea  junction,  at  the  intrusion  of  salt  water  under  fresh  water 
when  a lock  gate  is  opened,  in  the  ocean,  in  a collapsing  wake 
of  intermediate  density  on  an  interface  of  two  fluids  with 
different  densities  (in  other  words,  in  an  ocean  thermocline) 
and  finally  in  our  daily  lives  (thin  film  flow  on  an  inclined 
bed)  . 

In  hydraulics  this  is  called  the  lock  exchange  problem. 
Many  experiments  have  been  made  in  this  area.  A summary  can 
be  found  in  Turner  (1973)  . Benjamin  (1968)  showed  that  the 
front  must  have  a shape  of  head  behind  which  there  is  a turbulent 
region  and  an  abrupt  drop  to  a layer  of  uniform  depth.  K^rmin 
(1940)  showed  that  the  shape  of  the  nose  or  head  at  the  front 
is  60°  to  the  horizontal. 

Figure  6.1  shows  the  initial  geometry  and  its  velocity 
vector  plot.  A circular  cylinder  of  fluid  of  intermediate 


density  is  formed  by,  for  example,  the  propeller  of  a submerged 
vehicle  on  an  ocean  thermocline.  The  lines  show  the  location 
of  the  vertices:  171  points  altogether,  distributed  non- 

uni formly  over  the  first  quadrant,  with  the  higher  number 
density  near  the  thermocline  and  fewer  on  the  top.  This  is 
necessary  to  ensure  good  resolution  of  the  nose  geometry. 

The  vorticity  is  initially  zero.  Then  Equations  ( 2 S ) through 
(31)  are  applied  to  advance  the  calculation.  Due  to  the  lower 
fluid  density  over  the  cylinder  and  higher  underneath,  the 
buoyant  force  will  flatten  the  cylinder;  if  there  is  no  vorticity 
generated,  the  circular  cylinder  will  simply  be  flattened  into 
a thin  layer.  Due  to  the  vorticity  generated  by  the  buoyance, 
there  forms  an  advancing  nose  which  is  called  the  gravity 
current  or  weather  front  in  meteorology.  Notice  that  the  maximum 
velocity  at  t = 0 is  small.  As  it  developes,  Figure  6.2 
shows  at  t = .6  the  flattening  wake  and  its  velocity  distri- 
bution. The  velocity  has  grown  to  1.87  in  terms  of  the 
variables  defined  by  III.  2.  Figure  6.3  shows  the  well  defined 
nose  shape  at  t = 1.6.  The  nose  has  a slope  of  nearly  60° 
half  included  angle  as  predicted  by  Karman  (1940).  The  nose 
advancing  velocity  is  bounded  by 
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FIGURE  6.1b. 


Initial  Geometry  of  a Cylindrical 
Wake  at  Ocean  Thermocline 


FIGURE  6.3b 


Collapsing  Wake  on  an  Ocean  Thermocline 
at  t = 1.6 
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FIGURE  6.4a  Collapsing  Wake  on  Ocean  Thermocline 
at  T ■ 2.2 
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When  the  ratio  of  the  intruded  layer  depth  H to  the  overlying 
layer  depth  d,  H is  defined  by 


At  t = 1.6,  one  can  estimate  from  the  above  relations 
in  terms  of  the  presently-defined  non-dimensional  variables 
that  the  maximum  velocity  (the  velocity  at  the  nose)  is 
approximately  2 /tt  , which  agrees  with  what  is  shown  in  the 
velocity  vector  plots.  At  x = 2.2,  the  solution  is  shown  in 
Figure  6.4.  Notice  that  the  Kelvin-Helmhol tz  wave  develops 
on  the  lee  side  of  the  nose.  Although  this  is  a result  of 
the  numerical  noise  generated  by  the  finite  mode  approximation 
and  may  not  exist  in  reality,  it  certainly  is  a simple  mani- 
festation of  the  instability  of  the  flow  at  the  shear  inter- 
face. High  power  radar  probing  into  the  lee  of  weather  fronts 
did  find  a braided-like  structure.  In  this  connection,  we 
suggest  it  is  the  same  result  as  the  stationary  mountain  lee 
wave;  in  this  case,  the  mountain  is  replaced  by  the  moving  nose. 

I I I . 7 Buoyant  Wake  Near  an  Ocean  Thermocline 

In  this  section  we  turn  to  a simple  exten  ion  of  the 
subject  treated  in  section  III.  6;  -chat  is,  if  the  wake  is 
placed  at  a c_rtain  distance  below  a thermocline,  the  distance 
determines  the  upward  momentum  gained  through  the  rise  and 
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also  the  height  that  the  wake  overshoots  its  position 
of  neutral  equilibrium  - the  thermocline.  The  wake  is  assumed 
to  have  an  intermediate  density  between  the  fluids  above 
and  below  the  thermocline.  Due  to  the  negative  buoyance 
once  the  wake  overshoots  the  thermocline,  the  wake  will  be 
flattened  and  unless  it  falls  right  on  the  thermocline,  the 
vorticity  will  reverse  its  sign  alternatir.gly , 
while  spreading  itself  laterally.  In  the  absence  of  density 
stratification,  the  wake  will  simply  diffuse  and  grow  both 
in  vertical  and  horizontal  directions.  However,  when  the 
density  effect  dominates,  the  wake  will  be  squashed  and 
suppressed  so  that  the  residual  motion  consists  only  of  large 
scale  periodic  motion;  that  is,  the  internal  waves.  These 
waves  will  propagate  with  a typical  phase  velocity 
where  R'  is  the  wake  radius  when  it  overshoots  the  thermocline 
Since  R'  grows  with  t linearly,  it  is  expected  that  the  longer 
waves  generated  at  later  times  will  overpass  the  shorter  waves 
sent  out  earlier.  Wave  breaking  should  occur  and  a finite 
amplitude  wave  front  should  exist  as  the  case  in  III. 6. 

Figure  7.1  shows  tne  initial  geometry  of  a buoyant 
wake  underneath  a thermocline  at  ^ = 1.2.  Notice  that  the 
maximum  velocity  is  at  the  center  line  and  is  vertical  upwards 
At  t = . 6 , as  shown  in  Figure  7.2,  the  wake  ascended 
through  a vertical  distance  in  the  order  of  R = 1.  A well- 
defined  torus  is  formed  at  t = 1.38  shown  in  Figure  7.3. 


Since  it  has  overshot  the  thermocline,  the  negative  buoyancy 
force  is  in  action  to  flatten  the  top, and  vorticity  also  becomes 
smaller  on  the  top  portion.  One  can  estimate  the  vertical 
velocity  at  the  moment  using  the  arguments  presented  in 
Section  III. 2.  The  total  vortical  momentum  gained  through 
ascent  (a  distance  H = 1.5R)  is 

2 

A (mv)  = gAp  • A t ttP 
• ~ 2 

Assuming  m = pR  per  unit  length  in  the  z-direction,  and 

li 

substituting  At  — , we  get 


v • 

In  terms  of  the  length  R and  time  ^ we  qet  the 

I p 2 7T  R J 3 

non-dimensionalized  vertical  velocity 


V = 


which  agrees  with  the  value  appearing  in  Figure  7.3 
velocity  vector  plot. 

In  Figure  7.4,  the  outgoing  front  bears  some 
resemblance  to  that  in  Figure  6.3.  The  velocity  vector 
plot  shows  that  the  well-defined  vortex  is  skewed.  The 
horizontal  velocity  grows  at  the  horizontal  front.  One 
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FIGURE  7.1b  Buoyant  Wake  Near  an  Ocean  Thermocli 


FIGUkE  7.2a  Buoyant  Wake  Near  an  Ocean  Therraocline 
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FIGURE  7.2b  Buoyant  Wake  Near  an  Ocean  Thermociine 
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t h i now  ulocity  vacroa  plot 
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FIGURE  7.3b.  Buoyant  Wake  Near  an  Ocean 
Thermocline  at  x = 1.38 
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FIGURE  7.5b.  Buoyant  Wake  Near  an  Ocean 
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cannot  expect  the  solution  in  this  case  to  evolve  eventually 
into  the  one  shown  in  Figure  6.3,  because  during  the  rise 
the  cylinder  is  greatly  distorted  and,  in  the  absence  of 
viscosity,  the  wake  will  oscillate  indefinitely  on  the 
thermocline  at  the  Brunt-Vaisala  frequency.  The  horizontal 
velocity  at  the  advancing  front  may  again  be  shown  in  the 
right  order  of  magnitude,  a*  appears  in  Figure  III. 7. 4b. 

Figure  III. 7. 5 shows  the  final  plot  of  the  calculation  at 
t = 3 . 09 . The  velocity  vector  plot  indicates  there  is,  in  addition 
to  the  main  vortex,  a secondary  vortex  in  the  same  sense  of 
rotation.  The  solution  becomes  tortuous  and  further  calcu- 
lation seems  unrewarding.  One  can  always  reduce  the  time 
step  At,  which  was  determined  by  the  criterion  that  At  — ^ ^ 
~.03.  However,  the  calculation  will  be  costly  and  only 
the  smail  scale  motion  will  be  significantly  improved. 

HI-8  Saffman-Taylor  Instability 

Long,  narrow  convecting  cells,  that  is,  the  "salt  fingers", 
are  commonly  observed  when  hot  salty  water  is  poured  over  cold 
fresh  water.  A very  similar  phenomenon  occurs  at  the  interface 
of  two  superposed  viscous  fluids  when  they  are  forced  by 
gravity  and  an  imposed  pressure  gradient  through  a porous 
medium.  The  practical  examples,  in  addition  to  those  already 
mentioned  in  the  introduction,  are  oil-water  interface  in 
sand  or  in  shale  and  fresh  air-smoke  interface  in  a peat  moss 
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or  a granular  coal  bed  fire.  Saffman  and  Taylor  (1958) 
studied  the  finger-like  aLructure  in  a Hele-Shaw  cell  and 
found  chat  the  ratio  of  the  width  of  the  finger  to  the 
spacing  of  the  fingers  is  almost  always  equal  to  1/2. 

The  general  equation  [Equation  (21) ] is  composed  of 
two  diffusive  gradients.  A flew  system  in  this  context  is 
usually  called  the  doubly  diffusive  convection.  One  typical 
example  is  hot,  salty  water  overlying  cold,  fresh  water. 

Equation  (23)  determines  the  circulation  at  the 
discrete  vortex  (x^y.).  From  that,  the  velocity  field  is 
calculated  from  Equation  (14).  To  reduce  these  equations 
into  dimensionless  forms,  we  should  notice  that  the  flow  is 

characterized  by  two  quantities:  the  acceleration  2A£  ancj 

k P 

the  time  - , where  k is  the  permeability  and  n is  the  kinematic 

viscosity.  The  time  scale  is  derived  as  the  time  that  it 

takes  the  viscosity  to  diffuse  across  the  void  area  in  a 

porous  medium  which  is  represented  by  k.  From  these  two 

variables,  we  can  get  the  following  characteristic  dimensions: 

length  R,  k 2M  (!l)2 

time  T : — 

n 

circulation  ro=  i-  W2  ^ . 
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Libation  (23)  becomes 


( 

V 


u + i) An . 

O 1 


(33) 


U y. 

•/here  U is  the  dimensionless  variable  — T ami  u = — 
r.  K i R ' 

y _ 1 

'i  . Equation  (14)  is  reduced  to  Equations  (24)  and  (25) 

In  order  to  determine  the  maximum  allowable  time  step 
At,  it  is  necessary  to  find  the  order  of  magnitude  of  the 
terminal  velocity.  We  shall  attempt  to  estimate  this  quantity 
by  two  means. 

consider  that  the  finger  is  replaced  by  a sphere 

of  fluid  accelerated  under  the  effective  gravitational  force 

[Ap 


p and  decelerated  by  the  viscous  force  pV2u  — u Ute™inal 
times  the  surface  area  of  the  sphere.  Therefore 

- 4 aAp  k 

t 3 p n 

In  terms  of  the  characteristic  length  and  time,  we  have 
t 3 


Second,  if  we  assume  UQ  = 0 for  simplicity,  the 
vorticity  generated  is  that  due  to  the  terms  x ^ only. 

From  Figure  8*la,  the  vorticity  is  maximum  at  where  ^ 
is  the  maximum.  The  resultant  motion  will  be  to  lift  up  the 
center  and  push  down  the  external  edges.  At  some  later  time. 


114 


w V 1 1 - — ^ r 


Figure  8.1b  shows  the  growing  finger  where  the  vorticity 
of  the  opposite  sign  also  appears,  but  the  resultant  motion 
is  a further  acceleration  in  the  same  trend  as  in  Figure  8.1a. 
Assuming  the  final  stage  of  the  finger  structure  is  that 
depicted  in  Figure  8.1c,  we  can  estimate  the  velocity 
at  the  center  top  due  to  the  vortices  distributed  on  the  now 
vertical  interface  which  is  of  length  l. 


2 f*/2  _MiL.  . Jt  3M  ln  + II 

‘ VJo  ”*  P V !o 


2k  gAp 
hit  p 


i.n 


/i\ 

\h] 


if  h <<  ^ or  ut  = 9.2,  when  £ = lOh. 


From  these  two  estimations,  one  can  say  that  At  should 

be  approximately  i—  = .1  and  in  fact  we  set  At  = .01  to  be 

u t 

sure  of  a stable  time  integration. 

The  initial  disturbance  corresponds  to  that  shown  in 
Figure  8.2.  The  interface  is  perturbed  by  a gaussian 

displacement  at  the  center.  At  x = .22,  the  Figure  8-3 
shows  that  the  center  has  risen  while  the  edges  of  the  gaussian 
displacement  are  depressed.  Figure  8.4  shows  that  as  thf 
finger  grows  the  spacing  between  the  vortices  on  the  top 
becomes  large,  so  that  unless  a method  by  which  vortices  can  be 
added  to  this  region  is  implemented,  one  will  not  obtain  good 
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resolution  in  this  region.  Notice  that  in  Figure  8.5, 

at  t = .46,  one  vortex  has  been  added  into  the  center  region. 

Iliis  repacking  procedure  and  its  aspects  of  economy  are 
explained  as  follows. 

Tn  order  to  obtain  reasonably  good  resolutions  in 

describing  the  interface  boundary  without  being  committed  to 

using  a large  number  of  vortices  throughout  the  computation, 

it  is  essential  to  devise  a scheme  to  add  or  deplete  particles 

when  necessary.  At  initial  stages,  the  disturbance  slope  is 

small  and  the  separation  between  vortices  is  small.  Only  41 

particles  is  enough  to  yield  good  resolution.  As  the  fluid 

accelerates  toward  its  terminal  velocity  (defined  in  a manner 

similar  to  the  Stokes  terminal  velocity) , the  separation  grows 

so  that  additional  pat  tides  must  be  filled  in  wherever  the 

separation  is  too  large.  The  criteria  to  determine  the  repacking 

process  is  based  upon  the  separation  distance  between  two 

neighboring  particles.  Once  a preset  separation  is  exceeded, 

points  are  added  and  the  circulation  is  redistributed  among  the 

added  and  the  original  vortices  through  the  following  procedure. 

Assume  that  between  the  i^1  and  i + 1^^1  particles,  a 

particle  is  added.  The  new  particle  position  is  determined 

by  a linear  interpolation  between  ith  and  i+lth  vortices. 

The  circulation  on  the  new  particles  is  assumed  to  be  i (r.+T.  .) 

3 x l+l 

and  the  circulation  on  the  original  particles  is  reduced  to 

2 2 

j Fj  and  y ^i+1*  accordingly,  in  order  to  conserve  the  circulation. 
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FIGURE  8.3  Saf fman-Taylor  Instability  Study  at 
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FIGURE  8.5  Saf fman-Taylor  Instability  Study  at  T = 46 
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FIGURE  8.9.  Saffman-Ta 
at  i = .46 
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FIGURE  8.10.  Sa f fman-Taylor  Instability  Study 
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